Critical properties of the four-state Commutative Random Permutation Glassy Potts 

model in three and four dimensions 
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We investigate the critical properties of the four-state commutative random permutation glassy 
Potts model in three and four dimensions by means of Monte Carlo simulation and of a finite 
size scaling analysis. Thanks to the use of a field programmable gate array we have been able to 
thermalize a large number of samples of systems with large volume. This has allowed us to observe 
a spin-glass ordered phase in d = 4 and to study the critical properties of the transition. In d = 3, our 
results are consistent with the presence of a Kosterlitz-Thouless transition, but also with difi'erent 
scenarios: transient effects due to a value of the lower critical dimension slightly below 3 could be 
very important. 

PACS numbers: 75.10.Nr, 64.60.Fr, 05.10.Ln. 
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INTRODUCTION 



In the last years, spin-glass models without spin- 



inversion symmetry 



1.2.3.4.5.6.7.8.9 



have received a large 



amount of attention: probably the main reason for this 
big effort is that they are thought to describe structural 
glasses that in nature, as opposed to spin glasses, do 
not enjoy this symmetry. One of them, the Ferro-Potts- 
Glass (FPG),-i^ is a very direct generalization of the Ising 
Edwards- Anderson spin glass: the spins can take p dif- 
ferent values, and two neighboring spins contribute to 
the total energy a factor —Jij if they are in the same 
state and a factor +Jij if they are in different states. 
The bonds Jij are quenched random variables that can 
be distributed, for example, under a Gaussian or under a 
bimodal distribution. In the FPG, as we will discuss bet- 
ter in the following, the missing spin-inversion symmetry 
has the collateral effect of allowing the existence of a fer- 
romagnetic phase at low values of the temperature (this 
is why we define it Ferro-Potts-Glass) : because of this 
possible contamination the analysis of the glassy critical 
points of the model can potentially become very com- 
plex, and even lead to misleading conclusions. In fact, 
as we will discuss below, progress can be expected from 
the consideration of more refined models, where a gauge 
symmetry forbids the ferromagnetic phase. 

The FPG is a candidate for describing orientational 
glasses: a p-state spin models a quadrupole moment 
which can be directed in p (discrete) directions. How- 
ever, its main interest is maybe originated from some of 
the properties of its infinite-range version: for p > 4, 
for example, the mean field FPG undergoes a glass 
transition^ where the order parameter is discontinuous lii 
A number of different lattice modelSfi'^'^'^i^iiiSiii^ in other 



words, can be analyzed to clarify the finite-range behav- 
ior of systems showing the equilibrium properties typical 
of glasses: it is also important to remember that a num- 
ber of important connections have been foun d^'^'^'^ be- 
tween the mean-field dynamical equations of the model 
and the mode-coupling theory of the structural glass 
transitiou )^^'^^ that describes the evolution of the density 
correlations in a supercooled liquid above the dynamical 
transition temperature. 

Even if the mean-field results can be an important 
starting point, in a next step, since real systems have 
short-range interactions, it is important to study finite 
dimensional systems. Great part of the mainly numeri- 
cal effort has been focused on the p = 3 model in d = 3, 
to model a realistic quadrupolar glassiil The first nu- 
merical Studiesi^ii^i2ai21i2^i2^i2ii2^ found that the lower 
critical dimension di is close to 3. In a numerical study 
with a zero-temperature scaling approach, Banavar and 
Cieplalsi^ suggested that the FPG with Gaussian cou- 
plings has a di slightly greater than 3, while the FPG 
with bimodal couplings has a di slightly below 3 (but 
such a measurement had large intrinsic errors). A few 
months later Monte Carlo simulations^SiSi hinted that 
the transition seems to take place at a temperature com- 
patible with Tc = for both families of couplings, which 
suggested indeed that di = 3. Further simulations in the 
bimodal^^- and Gaussian^^ models were consistent with 
these results, although one could not exclude the possi- 
bility of Tc being small but larger than zero. A later study 
based on a high-temperature expansion^ii^ did not al- 
low to reach a final conclusion. Only recently we start to 
have clearer evidences about the situation: a large scale 
numerical study, based on a finite size scaling analysis of 
the correlation length-^, gives what looks like a reliable 
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evidence of a transition to a glass phase at finite Tc, mak- 
ing in this way a strong case for di being slightly below 
3 for the three-state FPG. 

Another interesting model that has been studied in de- 
tail is the p=10 model in d = 3, because of the intrinsic 
interest of the limit of a large number of states. Ol d^^i^^ 
and recent^finumerical simulations seem to suggest that 
there is no spin glass transition at finite temperature (but 
all the warnings about the dangers of ferromagnetic ef- 
fects at low T in this model stay in effect). This finding is 
in marked contrast with the predictions of mean field the- 
ory that indeed undergoes two transitions i^^'^^ new mod- 
els could be useful to understand better the connections 
among the mean field and the finite dimensional picture, 
and for example Potts-glass models with medium-range 
interactionsSii^ could be relevant at this effect. 

It has also been argued^^ (although some controversy 
exists^) that the choice of the coupling distribution 
might be relevant in removing the phase transition on 
the p= 10 model. The disease of the FPG that we have 
discussed before is the designated culprit: the lack of 
the spin inversion symmetry (which in Ising spin glasses 
is connected to a gauge symmetry that forbids a spon- 
taneous magnetization^) allows ferromagnetic ordering 
at low temperaturesi^i^ A partial relief to this problem 
can be obtained by using a distribution of couplings non- 
symmetric around zerOj-^v^^ but this choice does not re- 
cover the lost (important) gauge invariance. 

A different (and natural) definition of a frustrated 
Potts model containing quenched disorder, the Random 
Permutation Potts Glass (RPPG), has been introduced 
a few years agoi^ The key point of the RPPG (and of the 
similar model where only a set of possible couplings is 
allowed, the Commutative Random Permutation Potts 
Glass, CRPPG, where an additional symmetry is very 
useful to help checking thermalization, see III A[) is that 
it retains the gauge invariance which prevents Ising spin- 
glasses from entering ferromagnetic ordering at low tem- 
perature. The same paper® analyzed numerically the 
p — 4:, four-dimensional model (both in the RPPG and 
in the CRPPG versions) on lattices of volume V — A'^ 
and V — 5*. The two models were found to exhibit the 
same critical behavior, with a glassy phase characterized 
by a divergence of the overlap susceptibility. A prelimi- 
nary value of 7 was estimated from that divergence, and 
the critical temperature was obtained from the analysis 
of the Binder parameter: the critical behavior was found 
to be reached under a discontinuity, that was related to 
the one observed in the Random Energy Model. — It is 
also interesting to note that Carlucc:"^^ has discussed the 
relation connecting the (C)RPPG and the Chiral Potts 
model introduced by Nishimori and Stephen;^ which in 
mean field shows the same type of transition for p > 4i^!^ 
The authors of Ref. 6' also present a dynamical study of 
their models, and they observe clear aging effects. 

In this work we investigate, by means of Monte Carlo 
simulation and Finite-Size Scaling analysis, the critical 
properties of the three and of the four dimensional p = 4 



CRPPG. In d = 3, the finite-size behavior makes possible 
that the system undergoes a Kosterlitz-Thouless transi- 
tion, although a di barely lower than 3 is surely compat- 
ible with the significance of our numerical data. In d — A, 
we confirm the existence of the spin-glass transition re- 
ported in Ref. IH, but the use of a field programmable gate 
array (FPGA) computer (see the appendix and Ref. [33 ) 
allows us to obtain more accurate estimates of the criti- 
cal exponents, universal dimensionless quantities and non 
universal critical couplings of the model. 

The remaining part of this work is organized as fol- 
lows. In Section III Al we define the model and comment 
on its symmetries. We describe the relevant observables 
in Section IIIBI Section IIIII is devoted to a discussion 
of the numerical methods: the details of the simulations 
are described in Section [III Al and the Finite-Size Scaling 
method in Section fill Bl Further details about the com- 
putation are given in Section [ill CI while the problem of 
thermalization is addressed in Section [III DI The results 
for the d = 3 model are discussed in Section [TVl and those 
for d = 4 are discussed in Section [Vj We present our con- 
clusions in Section IVII In the Appendix we give details 
about the FPGA and about how they have actually been 
used. 



II. THE MODEL 

A. Model and symmetries 

We consider a system of spins {(Ti} defined on a d = 3 
(and d — A) dimensional simple cubic lattice of linear size 
L (volume V = L'^) and periodic boundary conditions. 
The Hamiltonian is: 

■H" = - XI ^'y.:n.,(a,) , (1) 
<»j> 

where the sum runs over all pairs of nearest neighboring 
sites. The spins can take the values {0,1,2,3}, and TTy 
are quenched permutations of {0, 1, 2, 3}, defined on the 
links of the latticei^i^ We define our quenched couplings 
(to implement the commutative model of Ref. [y) by 
extracting random permutations of (0, 1, 2,3) that com- 
mute with our "reference permutation" i?= (0, 1, 2, 3) — > 
(2,3,0,1). Only links from i to j such that ai — Uijlaj) 
give a non-zero contribution to the energy. The RPPG 
and CRPPG are deeply connected^ to the Chiral-Potts 
model analyzed by Nishimori and Stephen.— 

The symmetry with respect to the reference permuta- 
tion R helps in defining an order parameter q governed by 
a probability distribution symmetric under q —> —q (this 
turns out to be crucial for checking that the system has 
reached thermal equilibrium^). We define two copies of 

the system (two real replicas) {cp'}, {cp^} and we allow 
them to evolve independently at the same temperature 
and the same realization of quenched random couplings 
Ilij. The modified overlap between the two replicas at 
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site i is defined as 



III. NUMERICAL METHODS 



1 -f (1) (2) 

-1 if ^ ^ and erf ^ = (af ^ + 2) mod 2 



elsewhere . 



(2) 



B. Observables 



The main quantities that we will consider here are de- 
fined in terms of the Fourier transform of qi: 



qik) 



(3) 



The momentum space propagator is defined from the re- 
lation: 



G(fc) = y(g(fc)2) 



(4) 



In the thermodynamic limit and at the critical point, the 
propagator is expected to have poles at k — 0: 



G(fc) 



(5) 



where the correlation length ^ diverges at the critical 
point, and ^||^|| <C 1. We also define the non-connected 
susceptibility: 



X = G(0) . 



(6) 



On a finite lattice an extremely useful definition of the 
correlation length can be obtained from the discrete 
derivative of G{k). Using k = (27r/L)e^, where be- 
longs to the canonical Cartesian basis, one obtainsi^i^^ 



/ G(d)/G(fc) - 1~ 



1/2 



\^ 4sin2(7r/L) J 
We also compute and analyze the cumulant: 



C/4E 

We define the energy as: 
E 



(g(0)2)2 



(7) 



(8) 



(9) 



so that it lays in the [0, 1] interval. When we need 
to estimate the derivative with respect to /3 of an ob- 
servable O, we estimate it by measuring the connected 
correlation function {O H)c- Bias-corrected^S reweight- 
ing techniques^^i^i^^ allow us to use the numerical data 
taken at temperature T to compute expectation values 
at nearby temperature values T' , and to get in this way 
estimates that cover all the relevant part of the critical 
region. 



A. Simulations 

In the d — S model we have analyzed lattices of linear 
sizes L = 6, 8, 10 and 16. The critical behavior of the 
model (see Section livp has suggested to simulate a wide 
range of values of /3, ranging from 1.5 to 2.7. We have 
analyzed between 200 and 400 different samples of the 
smaller systems and around 1000 samples for L — 16. 

In d = 4, we have analyzed lattices of linear sizes L — 
8, 12, and 16, with (3 ranging from 1.385 to 1.5. The 
main computer effort has been accomplished around /3 = 
1.405 and (3 = 1.41, close to the critical point. At these 
temperatures, we have simulated 1000 samples for L = 8 
and 2000 samples for L> 8. For the other f3 values we 
have simulated between 200 and 400 samples. We have 
also analyzed 50 samples of the system deep into the low- 
temperature region, at /3 = 1.5. 



B. Finite size scaling 

We give here a few details about the finite size scaling 
approach that we have used for our analysis. When using 
the quotient methodM^^^^ one compares the mean value 
of an observable O, in two systems of sizes Li and L2, 
using the value (3 where the correlation length in units 
of the lattice sizes coincides for both systems. If, for the 
infinite volume system, {0){(3) oc |/3 — /3c|~^° , the basic 
equation of the quotient method is: 



io 



{0{(3,Li)) i{L2,m _L2 

£(il,/3) - ii 

it) ^'^^oLT 



(10) 



+ ...), 



where the dots stand for higher-order scaling corrections, 
V is the correlation length critical exponent, uj is the (uni- 
versal) first irrelevant critical exponent, and Aq is a non 
universal amplitude. 

Just below the lower critical dimension, at a distance 
e, the critical exponent \/v is expected to be of order 
e. This means that, for a limited range of lattice sizes, 
the slope of the ^/L curves at Tc grows very slowly (al- 
most logarithmically) with L. This could make life hard 
for a numerical study where one looks for a crossing of 
the ^/L curves, since the curves for the different lattice 
sizes would be basically parallel in the critical region. In 
other words, distinguishing a merging of the ^/L curves 
from a crossing becomes very hard. If one works pre- 
cisely at the lower critical dimension (i.e. e = 0), one 
may expect that one of two mutually excluding scenar- 
ios is realized. If Tc = 0, the curves for ^/L would not 
join (if plotted versus 1/T, the curves for lattices of size 
L and 2L should displace uniformly by a L-independent 
amount). On the other hand, if Tc > one would have 
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a Kosterlitz-Thouless picture, where the curves for ^/L 
merge for ah T < T^. It is clear that distinguishing a 
KosterUtz-Thouless scenario from e > but very smaU is 
numericaUy chaUenging. 

The most precise way of extracting the critical point /3c 
is to consider the crossing point of dimensionless quanti- 
ties such as ^/ L and U4. When comparing their values in 
two systems of size Li and L2, one finds that they take 
a common value at 



I3c + B 



{L2/L,y/--1 ' 



(11) 



The non universal amplitude B depends on the dimen- 
sionless quantity that one considers. 



C. Computational details 

In order to compute equilibrium expectation values we 
update the spins with a sequential Metropolis algorithm, 
we bring them to equilibrium and during the equilibrium 
dynamics we measure the interesting physical quantities. 
Thanks to our optimized FPGA based processor we have 
been able to run large scale simulations: for example 
thanks to strong thermalization tests we can be sure that 
we have thermalized systems of volume V = 16"^ and V = 
16^ at high /? values, already deep in the broken phase. 
We define an elementary Monte Carlo sweep (EMCS) 
as V sequential trial updates of lattice spin (considered 
in lexicographic order). To produce the needed pseudo- 
random numbers we use the Parisi-Rapuano shift register 
mcthod4i 
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5 X 10^ 


16 


2.4 
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600 


2 X lO'' 



TABLE I: For each lattice size of the d = 3 model, we show 
the simulated temperatures, number of samples, number of 
EMCS per sample and EMCS per measurement. 



The d = 3 small lattices, from L = 6 to 10, have 
been simulated at the cluster of the Instituto de Biocom- 
putacion y Fi'sica de Sistemas Complejos (BIFI). We have 
taken our measurements after every 40 EMCS. The total 
simulation time for this set of lattices has been equiva- 
lent of 0.2 years of a Pentium IV processor running at 
3.2 GHz. Our main effort in d = 3 has concerned the 
large, L = 16 lattice and has been simulated in a single 
FPGA (see lVII for details). The total simulation time cor- 
responds to almost 22 years of Pentium IV at 3.2 GHz. 
Table U shows the details of the computation. 

In the (i = 4 model, lattices with L = 8 and L—12 have 
been simulated at the BIFI Cluster. The total simulation 
time has been the equivalent to about 3 years of Pentium 
IV at 3.2 GHz. Again, the core of the simulation corre- 
sponds to lattice L = 16, and has been computed with 
the FPGA. The total simulation time has been about 
300 years-equivalent of Pentium IV. Measurements have 
been made every 5 x 10^ EMCS. The details of the com- 
putation are shown in Table HIl 



D. Thermalization tests 

This large computer effort has allowed us to thermalize 
in the broken phase lattices of volume including up to 
65536 spins (a large number). The thermalization issue is 
crucial in spin-glasses, and we have checked it by several 
independent tests. 

As a first tool we have used a logarithmic binning 
procedure. Let us say that during a Monte Carlo sim- 
ulation we have collected estimates for an observable 
quantity O at all integer times t in the interval [0, T). 
We divide these data in bins /„ = [T/2"+\T/2") for 
n = 0, 1, 2, 3, . . .. The usual disorder average of O, (O), is 
obtained (after assuming that all data are at equilibrium) 
by averaging all Monte Carlo data, i.e. the data over all 
bins. Information about thermalization can be obtained 
by averaging separately over samples the time series in 
the different bins. We get in this way the logarithmic 
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TABLE II: Same as Table H] for d = 4. 




FIG. 1: (Color online). Logarithmic data binning analysis 
(see text) of the non-connected susceptibility for the d — 4 
model, L = 16, /? = 1.41 and /? = 1.5. Notice that the large 
time region appears on the left in the figure. 



running disorder averages 0„ = (0)„. In usual logarith- 
mic data binning, if thermalization has been achieved, 
one expects that 0„ becomes rt-independent for small n 
(the last bins). We show this quantity (shifted by Oo 
for a better comparison with J„0, see below) in the case 
of the non-connected susceptibihty as a function of the 
logarithmic binning level n in Figure [TJ The data corre- 
spond to the four dimensional system of volume V = 16^, 
at two values of the temperature, one very close to the 
critical point and one in the low temperature phase: the 
errors are drawn with a thin line. 

An even better control of the convergence with time 
to the asymptotic result can be obtained by computing 
the difference of the thermal expectation value in bin n 
and the value in bin in each sample, and averaging 
this quantity over the disorder. In other words we de- 
fine 6nO = {0)n — {0)o- This way, one can obtain much 
smaller statistical uncertainty: we plot this quantity for 
the non-connected susceptibility in Figure [1] by drawing 
the errors with thick lines. 

For both f3 values of Figure[T]both indicators show that 
convergence has been reached. Errors in 6nX (thick error 
bars) are much smaller, but they still show that the last 
part of our samples has reached a steady state (even if 
the error is very small all the data of the last bin are at 
the level of one standard deviation from zero: also no- 
tice that the data for different data bins are correlated, 
that implies that correlated discrepancies have to be ex- 
pected) . We can claim that the data of the n = bin are 
surely well thermalized, and we use them for computing 
the equilibrium expectation values that we discuss in this 
note. 

We have also estimated the integrated autocorrelation 
time T for the observables that we have measured: we 
want to be sure that the total time length of our numer- 
ical simulation is far larger than r. 
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FIG. 2: (Color online). Correlation length in units of the 
linear size L as a function of /3 for d = ?> systems of different 
volumes. 



In (i = 3, for our larger system, L = 16, at fi = 2A (a high 
value of /3, deep inside the broken phase), we find that for 
the internal energy r = 5 x 10^ EMCS (and it turns out 
to be smaller for the other observables). This implies 
that our numerical simulation has been running for a 
time close to 12t. In d = 4, the length of the numerical 
simulation of the L — 16 system at f3 values close to the 
critical point turns out to be close to lOr. 

We have also used a further test of thermalization, by 
considering the data of the n = bin. We have done that 
by selecting a set of (3 values to use as starting points 
of the reweighting extrapolation. 3^. Figures [5] and [3] show 
an example of how data originated from different disor- 
der samples and independent numerical simulations yield 
consistent results. The choice of using different set of 
samples for different /3 values (the starting points of the 
different reweightings that appear in the figure as neigh- 
boring groups of points of the same type) does not opti- 
mize the quality of the final extrapolation of the data (in 
the full P interval that we consider), but gives a further 
check of both the quality of the thermalization and of 
the quality of the sample average. In our case the test is 
obviously successful. 

Even if these general thermalization checks are very 
useful, and they give strong hints that the system is ther- 
malized, the Z2 symmetry of the model (see scction lTl Ap . 
that has been introduced exactly with this goal in mind, 
is crucial to check thermalization. Let us repeat that 
the allowed couplings have been selected exactly such 
that the probability distribution of the modified overlap, 
P{q), has to be symmetric at equilibrium. We show in 
Figs.[3]and[n]P((7) for c? = 3 and d = A (computed by using 
the data of the n — bin, i.e. the last half of the data 
of the numerical simulation). These disorder averaged 
distributions show very clearly the expected symmetry. 

At last we have also studied the dynamics of different 
observables (for example of the modified overlap) in in- 
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FIG. 3: (Color online). The cumulant U4 as defined in Eq. |8] 
as a function of P for d = 3 systems of different volumes. 
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FIG. 4: (Color online). Distribution of the overlap in the 
d=3, I/=16 system, at several temperatures. 



dividual samples, and we show an example in Fig[51 We 
can observe a number of complete reversals of the global 
modified overlap, that gives us a new estimate of the time 
scale on which the system gets modified: this time scale 
is compatible with what we have estimated before. We 
stress again that the determination of this time scale is 
further evidence that we are indeed at thermal equilib- 
rium. 

We believe that this discussion clearly shows that it is 
safe to use for an equilibrium analysis the data from the 
n = bin (i.e. the last half of the simulation), since it is 
fully thermalized. 



IV. RESULTS FOR d = 3 MODEL 

We show in Fig. [2] the correlation length in units of L 
as a function of (3 for the three-dimensional model. In the 



FIG. 5: (Color online). Distribution of the overlap in the 
d = 4 model at low temperature (/3 = 1.44) for two different 
lattice sizes. 
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FIG. 6: (Color online). Evolution of the overlap of a repre- 
sentative sample of the d = 4 model, L = 16 system. Here 
P = 1.5. 



high-temperature regime the curves for different lattice 
sizes are well separated: for increasing (3 the different 
curves approach, and for values of (3 close to 2.3 they 
seem to have merged in a single curve. In the limits of our 
statistical accuracy, we do not see any sign of a splitting 
of the curves in the high-T phase (a crossing point at Tc 
and a splitting in both the low T and in the high-T phase 
is the usual signature of a usual phase transition): such a 
merging (without an eventual splitting) for increasing /3 
is what would happen in a Kosterlitz-Thouless transition 
(KT, see for example Ref. IH). 

The first (of the many) delicate issue about this poten- 
tial behavior concerns thermalization of the system: we 
have to be sure that we are not being mislead by the fact 
that we have not thermalized the larger lattice sizes (this 
could produce an effect hiding a crossing in the high-/? 
region). This is why we have studied, and discussed be- 
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fore, thermalization in detail: the thermalization checks 
described in Section IIIII make us confident that we have 
reached equihbrium for all the lattice sizes that we have 
considered. We should not forget that there are other 
possible issues that could hide from us, even in a very 
large scale simulation like the one discussed here, the 
asymptotic result: we could need for example a better 
statistical accuracy to discriminate a weak crossing, or 
we could need large lattices to see the crossing appear- 
ing, or we could need to go to higher /3 values. The issue 
of a very weak transition is a very delicate one, and re- 
liable statements must be phrased with great care. Here 
we claim that a KT scenario is a possible choice given 
the data that we have been able to measure in d~3, 

In a KT scenario the quantity ^/L is expected to re- 
main invariant in a finite low-temperature region adja- 
cent to the critical point. One way to be quantitative 
about that is to compute the crossing points /3^^'^^ for 
the dimensionless quantity U4, see Eq. [TTJ In Fig. [31 we 
plot the cumulant U4 for several lattice sizes. The curves 
for different lattice sizes cross close to /3 = 2.0 (look for 
example at the L = 8 and the L = 16 lattices) , at a temper- 
ature where the curves for ^/L on different lattice sizes 
did not yet merge (i.e. where the correlation length has 
the high-T behavior). The region of the crossing is quite 
narrow, so that is very implausible that the scaling cor- 
rections to U4 (usually larger than that of £,/L) will shift 
the crossings as much as to get them close to P — 2A. 
Therefore, under our numerical accuracy, we do observe 
that ^/L remains invariant in an interval of temperatures 
lower than that of the crossings of the cumulant. 

The features we have described are consistent with a 
transition of the KT type."*^ Nevertheless, as we have dis- 
cussed before, many possible effects could lead to difficult 
conclusions (for example the value of the lower critical di- 
mension to be slightly smaller than three). It is clear, in 
any case, that in d — S we are indeed sitting very close to 
the lower critical dimension. 




FIG. 7: (Color online). Correlation length in units of L as a 
function of /3 in the d = 4 model. 
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FIG. 8: (Color online). Zoom of the data of Fig. [7] close to 
the estimated critical point. 



V. RESULTS FOR THE d = 4 MODEL 

The authors of Ref. [fil, where the CRPPG model that 
we investigate here was proposed, found that the four- 
dimensional CRPPG undergoes a transition to a spin- 
glass phase at T«1.5 (by analyzing lattices of size L = 4 
and 5). 

In order to analyze the transition, we study here the 
scaling behavior of quantities as ^/L and C/4, that are ex- 
pected to be L-independent at the critical point. In Fig. [7] 
we plot the correlation length in units of the lattice size as 
a function of (3. The reweighting extrapolations of these 
quantities for pairs of lattices Li and L2 do intersect in 
the region around P — lAl. In order to be sure of the 
existence of the crossing we have thermalized lattices of 
linear size L = 8 and L = 16 deep in the low-temperature 
region: the normalized correlation length of the larger 
lattice is well above the one of the smaller lattice for f3 



values ranging from 1.44 to 1.5. 

In Fig. [5] we zoom the region closer to our putative 
crossing. In this region we have also thermalized lattice 
of linear size L = 12, and we include the L = 12 data in 
the figure and in our analysis. 

In Table IIIII we give the values of the crossing points 
^^1,^2 obtained by the crossing of the ^/L curves. Al- 
ready from Fig. [5] it is clear that the accuracy of the 
size-dependent estimates (3^^'^^ is not high enough to 
allow to estimate scaling corrections. This is since reach- 
ing thermal equilibrium for L > 16 was not in the scope 
of our numerical simulation (bound to run on a single 
FPGA chip), while lattices with linear size L < 8 would 
have probably been too small to show true asymptotic 
scaling corrections. 

Since the cumulant U4 scales like (,/L at the critical 
point, it might have played the same role than ^/L (by 
using Eq. Ilip . However, we find that it has much larger 
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scaling corrections than ^/L, and that these corrections 
shift the crossing points to higher temperatures, out of 
the range that we have analyzed (and where we believe 
the real asymptotic critical behavior can be observed). 
We have therefore not used U4 in our study of the critical 
point. Our results compare fairly with the ones obtained 
in Ref. by analyzing systems of linear sizes L — 4 and 
L = 5 (/3 must be renormalized since our Hamiltonian 
differs by a factor 2 from the one of Ref. [H) . 

To obtain the critical exponents we consider the oper- 
ators dpS, and x, whose associated exponents, see Eq.lTUl 
are a;^^^ — 1/ + 1 and x^. — j = v{2 — rf) . Taking the 
logarithm of the quotients of these expectation values at 
the crossing points of ^/L, we obtain the effective size- 
dependent exponents that we show in Table Uni We can 
summarize our best estimate for the d — A exponents as 
= 1.41(1), C/L = .47(2), V - 1.1(2), ry = -0.31(3) 
and 7 — 2.5(4): these error are statistical in nature and 
cannot, obviously, fully take care of the systematic ef- 
fects. 

As was happening in the determination of the value 
of the critical coupling, the estimated exponents lack the 
precision necessary for obtaining a reliable infinite vol- 
ume extrapolation. Ref. was quoting a value of 7 in 
the range between 1.3 and 1.5, obtained from the study 
of the overlap susceptibility in the warm phase of a lat- 
tice L = 8. Although our estimate is not very close to this 
value, it is clear that we are still dealing with lattice of 
intermediate size, and that a careful analysis of scaling 
corrections, that we hope will soon be possible, will prob- 
ably lead to reconcile these results. Our results should 
characterizes, if universality holds, the spin glass tran- 
sition to a Potts Glass, independently from the detailed 
model one selects. 

Finally, we also show in table IIIII the finite-size esti- 
mates of the universal quantity /L, i.e. ^/L evaluated 
at the critical coupling. 

VI. CONCLUSIONS 

We have presented a numerical study of the 4-state 
CRPPG in d = 4, and, for the first time, in (i = 3: we have 
used Monte Carlo simulations, reweighting techniques 
and a finite size scaling analysis. In d = 3 our evidence 
clearly shows that we are very close to the lower critical 





L2 


Pc,i/L 


C/L 






7 


8 


12 


1.41(1) 


0.47(2) 


1.1(1) 


-0.35(3) 


2.6(2) 


8 


16 


1.41(1) 


0.47(1) 


1.1(2) 


-0.33(2) 


2.5(4) 


12 


16 


1.41(1) 


0.46(2) 


1.0(4) 


-0.29(5) 


2.4(9) 



TABLE III: Our best estimates for the size dependent effec- 
tive critical coupling and for a number of universal quantities, 
as obtained from {L\,L2) pairs. 7 is obtained from the hy- 
perscaling relation 7 = i/(2 — 77). 



dimension, and suggests that a Kosterlitz-Thouless like 
behavior is possible, even if we could be dealing with a 
transient effect. In d = 4 we are able to collect a large 
number of thermalized samples for systems defined on 
large lattices, of linear size L = 16. Thanks to such a 
large scale numerical simulation we are able to qualify 
the spin-glass transition first found in Ref. 6, and we ob- 
tain size-dependent estimates of the critical coupling, of 
the critical exponents v and 77 and of the scale-invariant 
quantity 

In both cases, the use of a FPGA gives us the power 
needed to achieve thermalization, a target very ambitious 
for standard computers. We have been very careful in 
checking thermalization, and also thanks to the built-in 
symmetry of the CRPPG we have succeeded in this task. 
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THE FPGA DEVICE 

The problem of the glassy state, for example, is a typ- 
ical problem of very high complexity. A large (maybe 
infinite) number of time scales is involved, and numeri- 
cal simulations have to try to give hints about dynamics 
at very long times: very large correlation and thermal- 
ization times imply that, already on lattices of medium 
size, a huge computational effort is required. This is a 
typical situation where conventional computers could be 
not enough to do the job. 

The use of FPGA programmable chips for the sim- 
ulation of spin systems has been proposed several years 
ago^: conventional computers are not optimized towards 
the computational tasks relevant for our typical calcula- 
tion, and a FPGA can be programmed (at run time) in 
order to optimize the execution of the specific problem 
that one wants to solve. 

FPGA devices comes with numerous embedded and 
sizable memory blocks (RAM blocks), and thousands of 
configurable logic blocks with programmable intercon- 
nections. A configurable logic blocks can be programmed 
to perform complex logic operations and provide storage 
(fiip-flop registers) at the same time. 

A number of features that characterize our model are 
indeed optimal for being dealt with by a FPGA: we have 
discrete variables that can take a small number of values 
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(four for our p = A system) , and the interaction is local in 
physical space. The Metropolis algorithm and the ran- 
dom number generators discussed in Section IIII CI have 
been implemented in the FPGA in a very effective way. 

RAM blocks have a natural 2D (width x depth) grid 
structure. A 3D cubic matrix of bits can be obtained by 
stacking many of them, and access to all of them with 
the same memory address corresponds to addressing an 
entire plane in a 3D grid. We consider one such structure 
per each bit needed to represent fields (and interactions) 
defined on the sites of a simple cubic lattice. 

Locality of interactions (nearest neighbors) allows for 
a high grade of internal parallelism: in a checkerboard 
scheme, all black or all white sites of a lattice plane can 
be updated simultaneously (i. e. at the same clock cycle). 
Moreover, when simulating two real replicas and mixing 
black (white) sites of a system with white (black) ones 
of its replica, all sites in a plane can be processed in par- 
allel. Simultaneous local updates can then be performed 
by replicating small computation cells, each executing 
the few simple logical operations to compute local ener- 



gies, and including a 32 bit comparator for the Metropolis 
test. Precomputed transition probabilities (that allows to 
avoid lengthy computations of transcendental functions) 
are stored as several small look-up tables in configurable 
logic {distributed RAM), and addressed by the computed 
energy variations values (each look-up tables serves two 
distinct computation cells). The iterative processes in- 
volving 32 bit integer arithmetics for random number 
generators have also been parallelized, by cascading many 
32 bit integer adders and xors, and allowing for the gen- 
eration of hundreds of 32 bit random numbers per clock 
cycle. For further details, see Ref. [3^ 

We use the FPGA device Virtex 4/LX200, manufac- 
tured by Xilinx. Depending on lattice size and number 
of parallel updates (between 64 and 256) our designs run 
at clock speeds between 50 and 100 MHz. 

In Ref. [s^ its performances have been compared with 
the ones of a 3.2 GHz Pentium IV device: for the d — 
3 model the FPGA performs 1800 times faster than a 
Pentium, while this factor is 2300 in c? = 4. 
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